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ABSTRACT 

We develop an analytic model for the hierarchical correlation amplitudes 
S j,g( R ) = ?i,fl(- R )/S lfl ( R ) [where j = 3,4,5, and is the jth order 

connected moment of counts in spheres of radius R] of density peaks and dark 



ON 



matter halos in the quasi-linear regime. The statistical distributions of density 
peaks and dark matter halos within the initial density field (assumed Gaussian) 
are determined by the peak formalism of Bardeen et al. (1986) and by an 
extension of the Press-Schechter formalism, respectively. Modifications of these 
o distributions caused by gravitationally induced motions are treated using a 

q ■ spherical collapse model. We test our model against results for S^ g (R) and 

S4, g (R) from a variety of N-body simulations. The model works well for peaks 
even on scales where the second moment of mass (£ 2 ) is significantly greater 
than unity. The model also works successfully for halos that are identified earlier 
than the time when the moments are calculated. Because halos are spatially 
exclusive at the time of their identification, our model is only qualitatively 
correct for halos identified at the same time as the moments are calculated. For 
currently popular initial density spectra, the values of Sj )9 at R ~ 10 /i -1 Mpc 
are significantly smaller for both halos and peaks than those for the mass, 
unless the linear bias parameter b [defined by b 2 = £ 2 g( R ) / '£2 0^0 f° r l ar g e R ) 
is comparable to or less than unity. The Sj i9 depend only weakly on b for 
large b but increase rapidly with decreasing b at b ~ 1. Thus if galaxies are 
associated with peaks in the initial density field, or with dark halos formed at 
high redshifts, a measurement of Sj >g in the quasi-linear regime should determine 
whether galaxies are significantly biased relative to the mass. We use our 
model to interpret the observed high order correlation functions of galaxies and 
clusters. We find that if the values of Sj t9 for galaxies are as high as those given 
by the APM survey, then APM galaxies should not be significantly biased. 
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Subject headings: galaxies: clustering-galaxies: formation-cosmology: theory- 
dark matter 



1. INTRODUCTION 

A fundamental problem in cosmology is to understand how the spatial distribution 
of galaxies (and of galaxy clusters) is related to that of the underlying mass. In 
standard models of galaxy formation, it is assumed that galaxies form by the cooling and 
condensation of gas within dark matter halos (White & Rees 1978; White & Frenk 1991). 
Nongravitational processes which cannot be modeled reliably are likely to be critical in 
determining the properties of individual galaxies, yet they should have little effect on the 
formation and clustering of dark halos. As a result, the problem of galaxy biasing can be 
approached by first understanding how dark halos are distributed relative to the mass. 

Dark halos are highly nonlinear objects and their formation and clustering has usually 
been studied using N-body simulations (e.g. Frenk 1991; Gelb & Bertschinger 1994 and 
references therein). Such simulations are limited both in resolution and in dynamical range 
and can be difficult to interpret. Our understanding of their results could be substantially 
enhanced by simple physical models and the analytic approximations they provide. Mo & 
White (1996) have developed a model for the second order correlation functions of dark 
halos based on the Press- Schechter formalism (Press & Schechter 1974; hereafter PS). Here 
we extend their work to derive a model for the higher order correlations. We also work out 
a similar model for peaks in the initial density field. 

In the PS formalism, one defines dark matter halos at any given time as regions of the 
initial density field which just collapse at that time according to a spherical infall model. 
This formalism can be extended so that it predicts not only the mass function of dark halos, 
but also a wide range of other statistical properties of the hierarchical clustering process; 
comparisons with N-body simulation data show detailed agreement (e.g. Bond et al. 1991; 
Bower 1991; Kauffmann & White 1993; and particularly, Lacey & Cole, 1993, 1994). Mo 
& White (1996) showed that the PS formalism and its extensions can be used to construct 
a model for the spatial correlation of dark halos in hierarchical models. They used the 
standard PS formalism both to define dark halos from the initial density field, and to specify 
how their mean abundance within a large spherical region is modulated by the linear mass 
overdensity in that region. The gravitationally induced evolution of clustering was treated 
by assuming that each region evolves as if spherically symmetric. The model was found to 
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give an accurate description of the bias function, bh(R,5) = Sh(R,S)/S, where Sh(R,5) is 
the mean overdensity of halos within spheres which have radius R and mass overdensity 5. 
A simple extension provides a surprisingly accurate prediction for the halo-halo two-point 
correlation measured in N-body simulations (Mo & White 1996; Mo, Jing & White 1996). 
It is clearly interesting to see if such a model can also be constructed for the higher order 
correlations of dark halos. 

In the peak formalism, one assumes that galaxies form at those peaks of a suitably 
smoothed version of the initial density field which rise above some density threshold. 
Kaiser (1984) introduced this idea to show how the strong clustering of Abell clusters 
could result from the statistics of high peaks in a Gaussian initial field. This formalism 
was later developed extensively by Bardeen et al. (1986, hereafter BBKS). These authors 
showed that if galaxies can be associated with high peaks of the initial density field then 
they should be more strongly clustered than the mass, an effect usually called "galaxy 
biasing" . Unfortunately it is not known how well galaxies correspond to high peaks of 
the initial field, and it is unclear how to deal with the problem that the present-day 
clustering of peaks differs substantially from that in the initial (Lagrangian) space because 
of gravitationally induced motions. As a result, the predictions of the theory have only been 
checked quantitatively by direct N-body simulation of the motion of the material associated 
with density peaks (e.g. Davis et al. 1985; Frenk et al. 1988; Katz, Quinn & Gelb 1993). 
In this paper we will show that the idea underlying the model of Mo & White can also be 
used to construct a model for the correlation functions of peaks in physical space. We use 
the theory of Bardeen et al. to define peaks from the initial density field, and to specify 
how their mean abundance within a large spherical region is modulated by the linear mass 
overdensity in that region. The gravitationally induced evolution of clustering is then 
treated, as in the model of Mo & White, by assuming that each region evolves according to 
a spherical collapse model. 

We describe our model in §2, and present detailed tests of its predictions against 
N-body simulations in §3. Then in §4 we demonstrate how it can be used to interpret the 
observed high order correlation functions of galaxies and of clusters of galaxies, and to 
determine whether or not galaxies are biased relative to mass. Finally, in §5 we summarize 
our main results. 



2. THE MODEL 

To calculate the high order moments of peaks and halos (together called "galaxies"), we 
use the general formalism developed by Fry & Gaztanaga (1993). Consider the present-day 
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mass overdensity field 5(x) = [p(x) — p]/p, where p(x) is the local mass density and p the 
mean density. When smoothed in spherical window W(x; R) with characteristic radius R, 
5(x) gives rise to a smoothed density field: 

o(x;it>) = Jw(\x-y\;R)5(y)d 3 y. (1) 

For a top-hat window, <5(x; R) is just the volume average of <5(x) within a sphere of radius 
R. The statistical properties of <5(x; R) are described by the one-point distribution function 
of the density field. Now let us denote the overabundance of "galaxies" in the same window 
by 5 g (x;R). If 5 g (~x;R) is completely determined by 5(x;i?), then we can write S g as a 
function of 5, 5 g = F(S), which should not depend on x. (Note that this can hold at 
best approximately; see below.) To simplify notation, we will omit writing explicitly the 
smoothing radius R associated with 5 and 5 g . In general, we can expand F in Taylor series: 

OO L 

* 9 = F(5) = £ (2) 

where bk are constant. Fry & Gaztanaga (1993) have shown that if the volume averaged 
j-point mass correlation functions £j(R) have the hierarchical form: 

Z J (R) = S&-\R), (3) 

then the transformation given by equation (2) preserves the hierarchical structure in the 
limit £ 2 (-R) *C 1- In this case one can write 

t s *(R) = Sj£?(R), (4) 

and for j — 3, 4 and 5, which are relevant for our later discussion, one has 

S 3 , 9 = b' 1 (S 3 + 3c 2 ), (5a) 

S 4 , g = b~ 2 (S 4 + 12c 2l S 3 + 4c 3 + 12c|), (56) 

S 5 + 20c 2 ^ 4 + 15c 2 S 3 2 + (30c 3 + 120^)^3 + 5c 4 + 60c 3 c 2 + 60^] , (5c) 

where = bk/b and b — b\. Thus to obtain the high order moments Sj >g for halos and 
peaks in the quasilinear regime, we need to work out the coefficients bk in the bias relation 
(equation 2). We will do this in §2.2 and §2.3. 

It is important to emphasize that the bias relation given by equation (2) is at best 
approximate. This is because 5 g in a spherical region with mean mass overdensity 5 must 
depend not only on 5 but also on the internal structure of the region and on the external 
tides acting on it. As a result, there must be scatter in the "galaxy" counts among spheres 



S 5 , g = b~ 
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with the same 5 and R. Such scatter should also contribute to the high order moments. As 
discussed in Mo & White (1996), this contribution must be important on scales which are 
not much larger than the linear Lagrangian sizes of the "galaxies" . On larger scales it may 
or may not be important, and since it cannot be estimated within our formalism, in the 
present paper we will assume it to be negligible and test the resulting predictions against 
N-body data. 



2.1. Initial density field 

The initial overdensity field 5(x) = [p(x) — p]/p is assumed to be Gaussian and so 
to be described by a power spectrum P(/c)5z)(k — ki) = (5(k)5(ki)), where 5(k) is the 
Fourier transform of 5(x) and 5o(k) is the Dirac delta function. We smooth the field <5(x) 
by convolving it with a spherically symmetric window function W(x; R) having comoving 
characteristic radius R (measured in current units). The smoothed field can be written as 

£(x; R) = J W(k; i?)<5 k exp(ik • x)d 3 k, (6) 

where W(k; R) is the Fourier transform of the window function W(x; R). Following BBKS, 
we define the order / moment of the smoothed field by 

af(R) = J W 2 (k; R)P{k)k 2l d?k. (7) 

The order zero moment, which we denote by A 2 (R), is just the rms fluctuation of mass in 
the smoothing window. 

For a given window function the smoothed field <5(x; R) is Gaussian and so has the 
following one-point distribution function 

p(5; R)d5 = ^1/2 ex P 

Since both S and A(R) grow with time in the same manner in linear perturbation theory, 
it is convenient to use their values linearly extrapolated to the present time. These 
extrapolated quantities will still obey equation (8). In the following, we write our formulae 
in terms of these extrapolated quantities. We will also omit writing explicitly the smoothing 
radius R, but often use subscripts to distinguish A, and other quantities, at different 
smoothing lengths [e.g. A = A(R ), A 1 = A(Ri)]. As another convention, we will always 
label the properties of "galaxies" using a subscript "1". The subscript "0" is reserved for 
the properties of larger uncollapsed spherical regions. 



2A 2 (R) 



d5 



A(R) 



(8) 
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2.2. High order moments of dark halos 



Mo & White (1996) have described in considerable detail how to derive the bias 
relation, Sh = F(S), for dark halos (so subscript "h" denotes "halos") in the PS formalism. 
Here we adopt their notation and work out the coefficients bk in equation (2) that are 
needed in the calculations of S^h- 

We assume that dark halos are spherically symmetric, virialized clumps of dark matter. 
The mass M x of a halo is then related to the comoving Lagrangian radius Ri of the region 
from which it formed by 



Mi 



An 



pRl 



(9) 



where p is the mean density of the universe. In an Einstein-de Sitter universe (which we 
assume throughout this paper), a spherical perturbation of linear overdensity Si collapses 
at redshift z\ = Si/S c — 1, where the critical overdensity for collapse S c = 1.686. According 
to the PS formalism, the comoving number density of halos, expressed in current units, as 
a function of M x and Z\ is: 



n(Mi,Zi)dMi = - - 



1/2 



p Si dlnAi 
Mi ~K[d\nMi 



cxp 



si 



2A? 



dMi 

"m7 



(10) 



Notice that in this formalism a class of halos must be defined by specifying both their mass 
Mi (or equivalently Ri or A x ) and their redshift of identification z\ (or equivalently Si). 

To derive the bias relation, we need formulae which relate halo abundances to the 
density field on larger scales. Bower (1991) and Bond et al. (1991) extend the original PS 
formalism to show that the fraction of the mass in a region of Lagrangian radius i?o and 
linear overdensity S which at redshift z\ is contained in dark halos of mass Mi (where by 
definition M x < M ) is given by 



r/A 2 

/(AiA|AoA)^Mi 



Si - S 



(2tt)V2 (A?- A§)3/ 2 



cxp 



(Si - S ) 2 
2(A 2 -A 2 ) 



aWi 



dMi. 



Thus the average number of M x halos identified at redshift z\ in a spherical region with 
comoving radius R and overdensity S is 



AT{l\0)dM 1 = ^-f(Ai,Si\A ,So)^-dMi. 



(12) 



Following Mo & White (1996), we obtain the physical space overabundance of halos in 
spheres which at the desired redshift z have radius R and overdensity S, using a spherical 
model. In such a model, each spherical shell moves as a unit and different shells do not 
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cross until they collapse through the zero radius. Thus the mass interior to each mass shell 
is a constant, which gives 

R = (1 + 6) 1 ' 3 R. (13) 

Furthermore, since dark halos in our PS model are defined to be objects identified at some 
specific redshift, the mean abundance of equation (12) can be taken as referring to halos 
of mass Mi identified at redshift Z\ within spheres of radius R(Ro, S , z±) and overdensity 
8(8o, zi). Under these assumptions the average overdensity of dark halos in spheres with 
current radius R and current mass overdensity 5 can be obtained from equations (10) and 

where V = AttR 3 /3, R = R(l + o") 1 / 3 , and 5 is determined from 5 by the spherical collapse 
model, as described in the Appendix. When considered as a function of 8, 5h in equation 
(14) just gives the bias relation. We assume that R 3> R\ so that A\ — Aq in equation (11) 
can be replaced by A\. Assuming also 5 <C 1 and using equation (A4) in the Appendix we 
can expand Sh in the form of equation (2). It turns out that the first five coefficients (which 
are relevant in our discussion) are 

b = 0, (15a) 
bi = 1 + ^V^, (156) 

h = 2(1 + a 2 )^± + ( V -f) 2 (v\ - 3), (15c) 



Si Vox 

h = 6(a 2 + a 3 )±_l + 3(1 + 2a 2 ) {^f (v\ - 3) + {^f ±zM±l^ (i 5d ) 

v 2 — 1 fv \ 2 

6 4 = 24(a 3 + a ±)^i— + 12 [°2 + 2(a 2 + a 3 )] [jj {v\ - 3) 

+4(1 + 3a 2 ) (|) 2 ±lM±l + (|) 2 _ 10 ^ + 15), (15e) 

where z/i = 0*1/ Ai; a 2 , a 3 and a 4 are the coefficients in the expansion of 5 (8) (see equation 
A4). Inserting bk into equation (5) we can obtain S 3) h (skewness), S^h (kurtosis) and S 5t h 
for dark halos. As we can see from equation (15), for a given z±, the high order moments 
depend on the mass Mi (or Ai) of the halos. When halos (identified at a given redshift 
zi) with a range of masses are considered, bk in equation (5) should be replaced by the 
values obtained by averaging them over Mi with a weighting of n(M 1; z±). The high order 
moments depend also on the dynamical evolution of the underlying mass density field, as 
manifested by the dependence on a 2 , a 3 and a^. However, as shown by Bernardeau (1992), 
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the 5 -S relation (which determines a^) given by the spherical collapse model, and the high 
order moments (Sj) of the mass distribution, depend only weakly on cosmological model in 
the quasilinear regime. Therefore the results for Sj t9 given by equations (5) and (15) should 
not depend significantly on a particular choice of cosmological parameters. 

Before moving on to the next subsection, let us consider some asymptotic properties 
of the bk in equation (15) and of the high order moments Shj derived from them. When 
5\ ^> 1 and V\ is not large, i.e. for small halos identified at early times, equation (15) gives 
b\ pa 1 and bk ~ for k > 1. It then follows from equation (5) that Sj t h = Sj, meaning that 
such halos are not biased relative to mass. In contrast, when v\ ^> 1 and 5\ is not large, i.e. 
for big halos identified at low redshift, we have bk = b\ for k > 1, and Sj^ are determined 
completely by the statistical properties of the initial density field, independent of both Sj 
and dk- The numerical values of the first several moments are: 



If |z/i| <C 5i and 5\ is not large, i.e. for small halos identified at low redshift, then 
b\ pa 1 — 1/Si and bk ~ — k\(dk-i + dk)/Si for k > 2. In this case Sj t h may depend 
significantly on the dynamical evolution of the underlying mass density field. The skewness 
of such halos will be ~ [5i/ (5i — 1)]S^ — 6<5i(l + a 2 ) / (Si — l) 2 , which can be larger than 
S3. For halos with ui = 1, the skewness is S 3; h — S3 — which is substantially smaller 
than £3 unless 5\ (and so zi) is high. 



The argument in §2.2 can also be used to construct a model for the high order moments 
of density peaks. In the peak theory, we consider peaks in the initial density field after 
smoothing with a (spherical) window function with a given radius R\ (corresponding to a 
rms mass fluctuation of Ai), and examine the distribution of peaks with respect to the 
peak height v\ = S1/A1. According to BBKS, the comoving differential peak density is 



£>3,h — 3, S Ayh — 16, £5^ — 125. 



(16) 



2.3. High order moments of density peaks 



1 



(17) 



e 



where 




7 



(18) 



a 2 {R l )a Q {R i y 



with cr Q , <7i, 02 defined in equation (7), and 




(19) 



- 9- 



with 



/(*) = 



r 



3a; 



erf 



1/2 



X 



+ erf - 



2 



1/2 



31x 2 




e -^/8 + l?_ _ _ 




(20) 



(see equation A19 in BBKS). To derive the bias relation for peaks, we also need formulae 
which relate peak number density to the mass density field on larger scales. Let us consider 
a spherical top-hat region with Lagrangian radius Rq and linear mass overdensity So. The 
number density of peaks (with characteristic radius R\ < Rq) in such a region is modulated 
by the background field, and assuming R x <C Rq, it can be written as 



1 



e u p /2 G(<y p ,>ypis p )diSp, 



where 



1p = 



7 



v\ - ez/ 



(21) 



(22) 



<l_ e 2)l/2' ~P (l_ e 2)l/2> 

with i/q = 5o/A and e = (viv Q ) oc A /Ai. The last relation for e follows from the argument 
of Bower (1992) that (SiS ) oc A 2 , when R 3> R\. Under the same assumptions made for 
equation (14), the average overabundance of peaks S p (subscript "p" for "peaks") in spheres 
with current radius R and current mass overdensity S can be obtained from equations (17) 
and (21): 



s P m = 



1, 



(23) 



ra(i/i)V 

where Vq/V = (1 + 5), and So is determined from S by the spherical collapse model. As 
we have done for dark halos, we now assume Ro ^> Ri (so that Ai A ) and 5< 1, and 
expand S p in the form of equation (2). It follows that the first five coefficients are 



bo = 0, 



h = l + 



vj + 9i 
Si ' 



b 2 = 2(1 + a 2 



(24a) 
(246) 
(24c) 



h = 6(a 2 + a 3 )^±^ + 3(1 + 2a 2 ) (|) 2 (V 2 - 1 + 2 9l + ^ , 



93 



v x - 3(1 - gx)v x - 3^i + 6 U 2 + -J 



(24d) 
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2 2 / 

W = 24(a 3 + ^^jf 1 + 12 [«2 + 2(a 2 + a 3 )] K -1 + 2^ 



2^2 

^i 2 



+4(1 + 3a 2 ) 



MM 
5 J 5i 



^-3(1-^)^-3^ + 6 



. 93 
92 + -f 



^-(6-4^)^ + 3-12^ + 12 




(24c) 



where a 2 , a 3 , a 4 are, as before, the coefficients in the expansion of 5 (8) (see equation A4), 
and 



5% = 



(25) 



fc! G(7, 7^1) <9|/ fc 

Since G( , y,y) and its derivatives involve only single integrations (see equation 19), it is 
straightforward to calculate bk in equation (24) numerically. As one can see from equation 
(24), for a given peak scale R\ the high order moments of peaks given by equation (5) 
depend on the peak height v 1 (or 5i). When peaks with a range of heights are considered, 
bk in equation (5) should be replaced by the values obtained by averaging them over v 1 with 
a weighting of n(vi). It is interesting to note that bk in equation (24) would have the same 
forms as those in equation (15), if g± = —1 and gk = for k > 1. Since g% have finite values 
for any realistic power spectra, it is clear from equation (24) that in general the high order 
moments of peaks have different asymptotic values from those of halos discussed in §2.2. 



3. TEST BY N-BODY SIMULATIONS 
3.1. Simulations 

We now test our analytic theory by comparison with the results from a series of large 
cosmological N-body simulations of Einstein-de Sitter universes. We use results for four 
different spectra. The first two have CDM-like forms, where the transfer functions are 
given by equation (G3) in BBKS, with the shape parameter r = Vlh equal to 0.5 and 0.2, 
respectively. We normalize the initial power spectra by specifying <r 8 , the linear rms mass 
fluctuation in a spherical top-hat window of radius 8 /i _1 Mpc. These simulations were 
performed using a particle-particle/particle- mesh (P 3 M) code with 128 3 particles and a 
force resolution of about 0.2 /i _1 Mpc. The simulation box size is 256 /i _1 Mpc for V = 0.2 
and 300 /i _1 Mpc for T = 0.5. The two-point correlation functions and mass functions of 
dark halos in these simulations have been analysed by Mo, Jing & White (1996). 

The other two are power-law spectra with n = —1.5 and —0.5. These simulations were 
performed using the P 3 M code described by Efstathiou et al. (1988) and are very similar to 
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the simulations of that paper. However, they are substantially larger (N = 10 6 ) and have 
higher resolution (gravitational softening length equal to L/2500, where L is the side of 
the fundamental cube of the periodic simulation region). The initial power spectrum was 
normalized as described by Efstathiou et al. (1988) and "time" is measured by expansion 
factor a since the start of the simulation (a = 1 for the initial conditions). Jain, Mo & 
White (1996) have used these (and some other) simulations to test the similarity scalings in 
the mass correlation functions and power spectra. These simulations were also used by Mo 
& White (1996) to study the two-point statistics of dark halos. 

3.2. Tests for density peaks 

The peaks considered in this paper are defined as those above a certain threshold v s in 
the primordial density field smoothed with a Gaussian window exp(— r 2 /2r 2 ). The window 
width is taken to be r s = 0.54 h^Mpc so that it is relevant for galactic-sized objects. We 
follow the prescription of White et al. (1987) to select peaks in the numerical simulations 
(see Jing et al. 1994 for a detailed description of our algorithm). The algorithm gives an 
expectation number of peaks for each simulation particle. Since this number is always less 
than 1 (i.e. each particle carries less than one peak), we select peaks by randomly culling 
simulation particles with a selection probability for each particle equal to its expectation 
number. The correlation functions are calculated by peak counts in spheres regularly placed 
on a 32 3 grid. Since the simulations are periodic, there are no difficulties with spheres 
overlapping the boundary of the simulated region. 

The circles in Figure 1 show our simulation results for the skewness of peak count 
S^, P {R) as a function of the radius R of the spherical counting cell. Results are shown 
for peaks with heights above the values indicated in the panels. The solid curves show 
the predictions of equation (5a) with b(= b\) and b 2 given by equations (24b) and (24c), 
respectively. In our model calculations we have used the values of S 3 (R) (the skewness 
of the mass distribution) estimated directly from the simulations. In practice Ss(R) in 
the quasilinear regime can also be obtained from the initial density spectrum (Fry 1984, 
Juszkiewicz et al. 1993; Lucchin et al. 1994; Bernardeau 1994; Baugh et al. 1995). For 
comparison, we plot the skewness of the mass distribution in the simulations as the dashed 
curves. It is clear from Fig.l that the theoretical model works well over a wide range of 
scales. The thick ticks on the horizontal axis mark the value of R where the second moment 
of the mass distribution £ 2 (-^) = 1- O ur model can work well even for £ 2 (-^) > lj this is 
particularly the case when the peak-height threshold is not very high. Since very high peaks 
are preferentially located in high density regions where nonlinear evolution of the density 
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field can be strong, our model fails at small R for such peaks. 

In Figure 2 we compare the kurtosis S^ P (R) of the density peaks in the simulations 
(circles) with our model predictions (solid curves). Here again the skewness S 3 (R) and 
kurtosis S^R) of the mass distribution used in the model (equation 5b) are estimated 
directly from the simulations. For comparison, we show S 4 (R) in Fig. 2 as the dashed 
curves. Figure 2 shows that our model also works reasonably well for S^ p over a wide range 
of scales. 



3.3. Tests for dark halos 

To define dark halos in the simulations, we will use the standard 'friends-of-friends' 
(FOF) group finder with a linkage length equal to 20% of the mean interparticle distance 
(e.g. Davis et al. 1985). This algorithm is easy to implement and has been extensively 
tested against the PS mass function (see Lacey & Cole 1994 for a careful discussion). 
The two-point correlation functions of the halos found by this algorithm have been tested 
against analytical theories of the kind we analyse here by Mo & White (1996) and Mo et al. 
(1996). The statistics in this section are based on counts of halos or of individual particles 
within spheres. When evaluating such statistics for the simulations we count objects within 
spheres centered on each grid point of a regular 30 3 cubic mesh. 

Figure 3 shows the comparison of our model prediction for the skewness of halos S^h 
(solid curves) against results obtained from simulations (circles). Results are shown for halos 
in different mass ranges. In Figs. 3a-3c, halos are selected at an earlier epoch (when the 
expansion factor is a — ai) than the one when the skewness is calculated (at a — a 2 > ai). 
In general, halos identified at a\ will, by a 2 , have increased their mass by accretion or lost 
their identity by merging. However, galaxies which were forming at their centers at a\ 
may still remain distinct at a 2 . Thus the results shown in Figs. 3a-3c may be relevant to 
galaxies. In the simulations the position of each halo at the later epoch is assumed to be 
that of the particle which was closest to its center at a\. The model predictions are obtained 
from equations (5) and (15) with 5\ taking the value 5\ = (a 2 /ai)5 c , and with the skewness 
of mass distribution S 3 estimated directly from the simulations. In this case, the agreement 
between the analytic model and the simulation results is reasonably good for large R where 
the second moment of mass distribution £ 2 ~ 1- (The values of R where £ 2 — 1 are marked 
by the thick ticks on the horizontal axis.) The skewness of halos shown in the figure 
does not change significantly with halo mass, because the dynamical range covered by the 
simulation is still too small to allow us to select halos with masses small enough to see 
such a dependence. For comparison, the crosses in Fig.3a show the simulation result for 
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halos which contain only about 7 particles. We see that the skewness of such small halos 
can indeed be as high as that for the mass, as predicted by our model. Unfortunately such 
small halos are very poorly sampled by our current simulation. For small R, two effects may 
cause our model to fail. Both the nonlinear evolution of the mass density field, and the fact 
that halos are spatially exclusive at the time of their identification in the simulations, may 
change halo clustering properties on small scales. Halo exclusion effects reduce the variance 
in the halo count to significantly below the Poisson value and also affect higher moments 
of the counts (see Mo & White 1996 for more detailed discussion). For the same reason, 
our model is less successful for halos which are identified at the epoch when the skewness is 
calculated (Fig.3d); halo exclusion effects are clearly more important in this case. 

In Figure 4 we compare the kurtosis S^ h (R) of dark halos in the simulations (circles) 
with our model predictions (solid curves). Results are shown for the same models as in 
Fig. 3. As before the skewness S 3 (R) and kurtosis S^R) of the mass distribution, which are 
required by the model, are estimated directly from the simulations. The values of S^R) 
are shown in Fig.4 as the dashed curves. This figure shows that our model predictions for 
S^h agree reasonably well with simulation results on large scales when halos are selected 
earlier than the epoch when S^h is analysed. For the reasons discussed above, our model 
is less successful on small scales and for halos identified at the epoch when the kurtosis is 
analysed. 



4. IMPLICATIONS FOR GALAXIES AND CLUSTERS 

The last section shows that our model for the skewness and kurtosis of density peaks 
and dark halos ( "galaxies" ) works reasonably well. In this section we demonstrate how this 
model can be used to interpret the observed high order correlation functions of real galaxies 
and clusters of galaxies. The assumption we make is, of course, that these objects are 
associated with initial density peaks or with the dark halos present at some given redshift. 

In Figure 5 we show predictions for the high order moments Sj t9 (j = 3, 4 and 5) of 
"galaxies" as a function of the linear bias parameter b = bi [defined by equations (15b) 
and (24b) for halos and peaks, respectively]. Each curve corresponds to a particular choice 
of Si, as parameterized by the value of zi — Si/S c — 1 given in the figure caption. As an 
example, we show results for a CDM-like spectrum with V = 0.2 and <t 8 = 1 and for a 
radius R = 10 ft, _1 Mpc. The spectrum chosen here is consistent with that given by the 
angular correlation functions of galaxies in the APM survey (Efstathiou, Sutherland & 
Maddox 1990; see also Maddox et al. 1996 for a recent discussion). The main features in 
Sj t9 do not change significantly if we change the value of T from 0.2 to 0.5. The choice of 
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R is based on the fact that the mass density field in the universe is still in the quasilinear 
regime on this scales and that high order moments of galaxies are difficult to measure on 
much larger scales. For halos with fixed z x and b, Sj^R) depend only on Sj(R) which, in 
turn, depend only on the effective power index at R: n c g(R) = —3 — 2m(A)/mi?. For 
peaks, however, Sj jP (R) depend, in addition, also on the shape of the spectrum on the peak 
scales Ri [through the dependence on 7 defined in equation 18]. For a given shape of power 
spectrum, the dependence of Sj iP on a 8 is weak, once b and 5i are fixed. The values of Sj t9 
in the quasilinear regime depend only weakly on cosmology, as pointed out in §2.2. 

From Fig. 5 we see that the Sj t9 for "galaxies" are systematically smaller than those for 
the mass, unless the linear bias parameter b is comparable to or less than unity. The values 
of Sj t9 are more or less constant for large b, but decrease rapidly with increasing b for b £ 1. 
For a given b, Sj t9 are higher for dark halos that are identified at an earlier epoch and for 
density peaks with higher 5i. The values of Sj j9 are the lowest for "galaxies" with z\ — 
and b ~ 1. These results have interesting implications for the observed high order moments 
of galaxies and clusters of galaxies. 

Let us start with clusters of galaxies. Clusters of galaxies are found to have much 
larger two-point correlation amplitude than galaxies (see e.g. Bahcall & West 1992; Peacock 
& West 1992; Dalton et al. 1992; Nichol et al. 1992). If the two-point correlation function 
of galaxies is neither significantly smaller than that of the mass nor larger than that 
of the mass by a factor exceeding three, the value of the linear bias factor for clusters 
should lie somewhere between 2 and 5. If we take clusters to be virialized halos identified 
at present time, Figure 5 shows that the skewness (Ss iC ) and kurtosis (S^) of clusters 
should be about 2 and 7, respectively. This is consistent with current observations. Based 
on the three-point correlation functions of Abell clusters, Jing & Zhang (1989) found 
Q ~ 0.6 which corresponds to S^ c ~ 2 (see also Plionis & Valdarnini 1994 and Cappi & 
Maurogordato 1995 who did count-in-cell analysis for Abell clusters and confirmed that 
Ss tC ~ 2). Recently, Gaztanaga, Croft & Dalton (1995) obtained S^ c ~ 2 and S^ c ~ 8 
for clusters in the APM survey. These values of skewness and kurtosis appear to be much 
smaller than the corresponding values for the mass distribution obtained from a power 
spectrum which has the shape expected given the angular correlation functions of APM 
galaxies (see e.g. Gaztanaga et al. 1995). From our model we see that such low values are 
a result of clusters being high mass halos identified at low redshift. 

For optical galaxies, the current best estimates of the high order moments are those 
of Gaztanaga (1994) based on APM survey. The values he got are Ss t9 (R) = 3.16 ± 0.14, 
S A>g {R) = 20.6 ± 2.6 and S 5 , g (R) = 180 ± 34 for R ~ 10 h^Mpc. These results are plotted 
in Fig. 5 as circles with errorbars. As noticed by Gaztanaga & Frieman (1994), the observed 
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values for APM galaxies appear to be close to those for the mass (the values predicted 
by quasilinear theory are indicated by the horizontal solid lines in Fig. 5). Comparing our 
predictions with the observational results we infer that galaxies in the APM survey should 
not be strongly biased relative to the mass. Namely, the linear bias parameter for APM 
galaxies should not be significantly larger than unity. Mild antibias (i.e. b ^ 1) is possible 
if most APM galaxies are associated with galactic-sized halos which form rather late. We 
note, however, that Gaztanaga (1992) found significantly smaller values Ss, g (R) ~ 2 and 
S4, g (R) « 5 from the fully 3-dimensional CfA and SSRS surveys. Gaztanaga (1994) argues 
that these small values reflect the fact that the volumes of the local surveys are too small to 
be fair. Redshift distortion present in these surveys may also complicate the determination 
of the high-order correlations in real space. Nevertheless, if these lower values turn out 
to be correct (for example if there is some problem in deriving S^g and S^ g from the 
2-dimensional APM data) then our analysis would suggest that the APM galaxies could 
be substantially biased relative to the mass. Future redshift samples from large digital sky 
surveys will certainly help to resolve the problem. 

The skewness and kurtosis of spiral galaxies and IRAS galaxies appear to be much 
lower than those for APM galaxies, with S3 ~ 2 and £4 ~ 10 (e.g. Jing et al. 1991; Meiksin 
et al. 1992; Bouchet et al. 1993). Since these galaxies have weaker two-point correlations 
than APM galaxies (and therefore lower b values), the observational results seem to require 
these galaxies to be associated with halos identified at late time, or with peaks at low 
overdensity 5 1 . 

5. CONCLUSIONS 

In this paper, we have developed an analytic model for the high-order moments of 
the distribution of density peaks and dark halos in the quasi-linear regime. Such a model 
allows the high-order correlation functions of density peaks and dark halos to be calculated 
analytically for any given initial (Gaussian) density spectrum. Tests against results from 
a variety of N-body simulations have shown that our model works successfully for density 
peaks and for halos identified at an earlier epoch than the time when the moments are 
calculated. Our model is only qualitatively correct for halos identified at the same time 
as the moments are calculated, because halos are spatially exclusive at the time of their 
identification. We have found that the skewness (S^ g ), kurtosis (S^g) and £5^ for both 
halos and peaks decrease rapidly with the linear bias parameter b (of these objects) for 
b ~ 1. Thus if galaxies are associated with peaks in the initial density field, or with dark 
halos formed at different redshifts, a measurement of Sj j9 (j = 3, 4, 5) of the galaxy 
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distribution in the quasilinear regime should allow us to determine whether or not galaxies 
are significantly biased relative to the mass. We have used our model to interpret the 
observed high order correlation functions of galaxies and clusters. We have found that if 
the values of Sj t9 for galaxies are indeed as high as those given by the APM survey, then 
APM galaxies should not be significantly biased. 

There is, however, a significant uncertainty in the comparison between our model 
predictions and observed galaxy distribution. At any given time massive dark halos may 
contain more than one galaxy and galactic-sized peaks in the initial density field may merge 
with each other to form a single galaxy. Thus the observed galaxies may not correspond 
uniquely to the centers of the halos present at any single epoch or to galactic-sized peaks in 
the initial density field. As a result, it is not straightforward to apply our results directly 
to galaxies. However, if more detailed modeling allows a prediction of how galaxies form 
in dark halos (e.g. White & Frenk 1991; Kauffmann, Nusser & Steinmetz 1996), our 
results can readily be extended to study the high-order moments of galaxy distribution. 
Such a study will also help us to assess the importance of nongravitational effects in the 
measurements we are suggesting here. 
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A. The spherical collapse model 

For a spherical perturbation in an Einstein-de Sitter universe, the physical radius R of 
a mass shell which had initial Lagrangian radius R and mean linear overdensity S is given 
for 5 > by (see Peebles 1980) 

R(Rq,5 q ,z) _ 3 1 - cos 9 

Ro ~ 10 \S \ ' 1 j 

1 _ 3x6 2 / 3 (fl-sinfl) 2 / 3 
1 + z ~ 20 \6 \ ' ^ ' 

For Sq < 0, we just replace (1 — cos#) in equation (Al) by (cosh# — 1) and (6 — sin#) in 
equation (A2) by (sinh# — 6). Without loss of generality, we assume z = at the time when 
the moments of halos and peaks are examined. Then S depends only on the present mass 
overdensity S = (R /R) 3 — 1. For \S\ 1, we can expand S (S) in power series: 

oo 

So = Y, a kS\ (A3) 

k=0 

where the first five coefficients (which are used in our model) are 

17 341 55805 

«o = 0; a 1 = l; o 2 = --; a 3 = — ; a 4 = (AA) 

(see Bernardeau 1992). As shown by Bernardeau, the So-S relation depends only very 
weakly on cosmological model in the quasilinear regime. 
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Fig. la. — The skewness of density peaks with different heights v predicted by our model 
(solid curves) compared with that derived from N-body simulations (circles). The dashed 
curves show the skewness of the mass density distribution in the simulation. Results are 
shown for the standard cold dark matter model with (fi, T, er 8 ) = (1,0.5,0.62). The thick 
ticks on the horizontal axis show the values of R where £ 2 (-^) = 1- 
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Fig. lb. — The same as Fig. la for a model with (Q, T, a 8 ) = (1, 0.5, 1.24). 
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Fig. lc. — The same as Fig. la for a model with (Q, T, a 8 ) = (1, 0.2, 0.5). 
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Fig. Id. — The same as Fig. la for a model with (Q, T, a 8 ) = (1, 0.2, 1). 
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Fig. 2a. — The kurtosis of density peaks with different heights v predicted by our model 
(solid curves) compared with that derived from N-body simulations (circles). The dashed 
curves show the kurtosis of the mass density distribution in the simulation. Results are 
shown for the standard cold dark matter model with (fi, T, a 8 ) = (1,0.5,0.62). The thick 
ticks on the horizontal axis show the values of R where £ 2 (-^) = 1- 
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Fig. 2b.— The same as Fig. 2a for a model with (fi, T, a 8 ) = (1, 0.5, 1.24). 
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Fig. 2c— The same as Fig. 2a for a model with (Q, T, a 8 ) = (1, 0.2, 0.5). 
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Fig. 2d.— The same as Fig. 2a for a model with (Q, T, a 8 ) = (1, 0.2, 1). 
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Fig. 3a. — The skewness of dark halos with different masses (indicated by N p , the number 
of particles contained in the halos) predicted by our model (solid curves) compared with 
that derived from N-body simulations (circles). The dashed curves show the skewness of 
the mass density distribution in the simulation. Results are shown for scale-free model with 
n = —1.5. Halos are selected at an earlier epoch (when the expansion factor a has the lower 
value indicated in the first panel) than when the skewness is calculated (at the epoch with 
the higher a). The thick ticks on the horizontal axis show the values of R where £ 2 (-^) = 1- 
For comparison, we also plot (as crosses) the simulation result for small halos to show the 
increase of skewness with decreasing halo mass for such halos. 
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Fig. 3b. — The same as Fig. 3a for halos selected at another epoch. 
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Fig. 3c. — The same as Fig. 3a for a model with n = 



-0.5. 
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Fig. 3d. — The same as Fig. 3a for halos selected at the same epoch as the one when skewness 
is calculated. 
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Fig. 4a. — The kurtosis of dark halos with different masses (indicated by N p , the number 
of particles contained in the halos) predicted by our model (solid curves) compared with 
that derived from N-body simulations (circles). The dashed curves show the kurtosis of the 
mass density distribution in the simulation. Results are shown for scale-free model with 
n — — 1.5. Halos are selected at an earlier epoch (when the expansion factor a has the lower 
value indicated in the first panel) than when the kurtosis is calculated (at the epoch with 
the higher a). The thick ticks on the horizontal axis show the values of R where £ 2 (-^) = 1- 
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Fig. 4b. — The same as Fig. 4a for halos selected at another epoch. 
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Fig. 4c. — The same as Fig. 4a for a model with n = 



-0.5. 
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Fig. 4d. — The same as Fig. 4a for halos selected at the same epoch as the one when kurtosis 
is calculated. 



-35 - 




Fig. 5. — Model predictions for the high order moments Sj, g (j = 3, 4 and 5) of halos (left 
panels) and peaks (right panels) at a radius r = 10 h~ 1 Mpc as a function of the linear bias 
parameter b = b\. Each curve shows result for a given 5i (the linear overdensity of halos 
and peaks). Results are shown for z\ = (#i/1.68 — 1) = 0, 0.5, 1, 2 and 4 (curves from 
bottom up). The horizontal lines show the values of Sj for the mass density field calculated 
from quasilinear theory, whereas the data points (plotted arbitrarily at b = 2) show the 
observational results for APM galaxies (Gaztanaga 1994). 



